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ABSTRACT 

To perform nuclear reactor simulations in a more realistic manner, the coupling scheme between neutronics 
and thermal-hydraulics was implemented in the HNET program for both steady-state and transient conditions. For 
simplicity, efficiency, and robustness, the matrix-free Newton/Krylov (MFNK) method was applied to the steady- 
state coupling calculation. In addition, the optimal perturbation size was adopted to further improve the convergence 
behavior of the MFNK. For the transient coupling simulation, the operator splitting method with a staggered time 
mesh was utilized to balance the computational cost and accuracy. Finally, VERA Problem 6 with power and boron 
perturbation and the NEACRP transient benchmark were simulated for analysis. The numerical results show that 
the MFNK method can outperform Picard iteration in terms of both efficiency and robustness for a wide range of 
problems. Furthermore, the reasonable agreement between the simulation results and the reference results for the 
NEACRP transient benchmark verifies the capability of predicting the behavior of the nuclear reactor. 
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1 Introduction 


With the growing computing power and requirement of safety analysis [1], there has been interest in the high- 
fidelity simulation of the nuclear reactor core [2—4]. The operational phenomena in the reactor to be simulated are 
very complex owing to the close interaction among various physical fields, particularly neutronics and thermal 
hydraulics [5, 6]. More specifically, the power density from the neutron flux significantly affects the thermal 
properties, which in turn has a powerful influence on neutronics through cross sections. Hence, to achieve a realistic 
prediction of reactor behavior under steady-state and transient conditions, the strong bond between neutronics and 
thermal-hydraulics must be considered. 

Because the steady-state solution is generally regarded as the initial condition for transients [7], the first priority 
is to address the steady-state coupling of neutronics and thermal hydraulics. For steady-state coupling calculations, 
the Picard iteration has been studied and used extensively in high-fidelity simulation codes. However, convergence 
may be difficult and inefficient when the nonlinearity of the problem increases significantly [8]. Hence, some 
approaches have emerged to address the limitations of Picard iteration, such as relaxation, partially converging 
physics [8], and Anderson acceleration [9]. Alternatively, the fully implicit Jacobian-free Newton-Krylov (JFNK) 
method [10], with fast convergence and good robustness, has come into focus. However, most high-fidelity 
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neutronics codes do not form the residuals required to implement JFNK. Therefore, the application of JFNK requires 
significant modifications to the algorithms of these well-developed codes. For simple implementation and good 
convergence performance, the MFNK method was proposed [11]. In contrast to the JFNK method, MFNK 
constructs residuals through the solutions of individual subsystems within the Picard iteration framework. Thus, 
MENK can possess the advantage of implicit coupling while avoiding the reformulation of the coupled subsystems 
into a larger nonlinear system. 

After the steady-state coupling calculation, the transient coupling simulation was emphasized. Although the 
implicit coupling scheme mentioned above has been applied to the nodal method for transient coupling calculation 
[12], most transient coupled systems still adopt the operator splitting (OS) method to compromise between 
computational burden and accuracy. Numerical instability may occur because of the explicit exchange of coupled 
information. Some improvements to the OS method have been developed, and the two main strategies are using 
staggered time mesh and higher-order time discretization of the nonlinear terms [12]. 

The high-fidelity Neutron Transport Program (HNET) is a 2D/1D transport code based on multilevel 
acceleration theory [13, 14]. The accuracy and efficiency of HNET have been verified against various steady-state 
and transient tests [15, 16], but only for problems without thermal-hydraulic feedback. To consider the interaction 
between neutronics and thermal hydraulics, a simplified internal thermal-hydraulic module was developed in HNET 
to provide thermal-hydraulic feedback. The MFNK method was adopted for the steady-state coupling calculation. 
Moreover, the OS method with a staggered time mesh was applied to the transient simulation. 

The remainder of this paper is organized as follows. In Sect. 2, the MFNK method and the estimation of the 
associated optimal perturbation size are introduced. Section 3 describes the implementation of the coupling 
neutronics with thermal hydraulics. An assessment of the MFNK method for the VERA core physics benchmark 
problem 6 [17] and the numerical results for the NEACRP transient benchmark [18] are presented in Section 4. 


Finally, conclusions are drawn in Sect. 5. 


2 Matrix Free Newton/Krylov method 


2.1 Derivation of MFNK method 

The MFNK method is based on the block Gauss-Seidel (BGS)-style Picard iteration. For the nuclear reactor 
nonlinear system considering neutronics and thermal hydraulics, the BGS-style Picard iteration alternatively solving 
each physical field can be written as in Eq. (1) and Eq. (2). 


xk! =F, (y ) (1) 


y =F, (x) 2) 
where vector x denotes the fuel temperature, coolant temperature, and coolant density in the thermal-hydraulic 
model, and y denotes the power density derived from the neutron fluxes. Fry and Fy represent the thermal-hydraulic 
solver and neutronics solver, respectively. Superscript k represents the iteration number. Because the coupled 
problem is implicit, the overline in Eq. (1) and Eq. (2) indicates that the corresponding variables should be 
transformed further. In this sense, Fru and Fy should be regarded as abstract solvers involving the total required 
operations, rather than only an actual solver. For brevity, the tilde is omitted from the expressions in the subsequent 
analysis. 

When the coupled problem converges, the solutions to Eq. (1) and Eq. (2) must be consistent, which implies 
that the nonlinear system can be expressed in terms of a residual form, as in Eq. (3). 


F(y') = y' -Fy (Fy (y*)) =0 (3) 


Subsequently, Newton method was adopted to resolve the nonlinear equation. By applying Taylor’s expansion 
to the right-hand side of Eq. (3) around y*, and ignoring the higher-order terms yields Eq. (4). 


F(y' + dy) = F(y")+ F'(y" dy =0 (4) 
Then, Eq. (4) can be rearranged as shown in Eq. (5), which is referred to as Newton iteration. 


J(y' dy =-F(y") (5) 


where J(y*) represents the Jacobian matrix. 

The linear system Eq. (5) can be solved using the Krylov subspace method, where the Jacobian matrix is only 
utilized for the action of the matrix-vector product. In this scheme, the Jacobian matrix-vector product can be 
approximated by the finite difference approximation to avoid assembling the Jacobian matrix explicitly, as shown 
in Eq. (6). 
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where e is the perturbation size and v is the basis vector of the Krylov iteration. 
Upon the dy, the power density is updated: 
k+ k 
yu =y +oy (7) 


The MFNK iteration is repeated until the convergence criteria are met. Subsequently, thermal-hydraulic results 
are obtained based on the converged value y. From the above analysis, it can be found that only the solution and not 
the residual of individual subsystems are necessary to implement MFNK. 

2.2 Optimal perturbation size 

The finite-difference approximation implemented in Eq. (6) results in truncation and rounding errors, as well 
as convenience. The truncation error is proportional to the perturbation size, whereas the rounding error is the 
opposite. Therefore, a balance should be found between the truncation error and the rounding error through the 
selection of the perturbation size. 

Xu [11] proposed a method to estimate the optimal perturbation size. Assuming that y, v, and € are given at 
the machine precision, the error ea introduced by the finite-difference approximation can be arranged into Eq. (8). 
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where y+ev is the finite precision representation of y+év and F is the function F evaluated for finite 


precision arithmetic. The terms on the right-hand side of Eq. (8) are referred to as truncation and round-off errors, 
respectively. 

After quantifying the error caused by the finite-difference approximation [11], the error upper bound can be 
obtained. Subsequently, the optimal perturbation size can be estimated by minimizing the upper error bound, as 
shown in Eq. (9). 
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where Emach is the machine precision, y represents the Lipschitz constant, 7 is the relative error in computing 


Fy(Frx(y)), and I denotes the L-2 norm in this study, except where otherwise noted. 


Note that the numeric values of |F (y)| and y are unavailable during the MFNK scheme process. To 


overcome this difficulty, the following approximation of Eq. (10) and Eq. (11) were adopted. 
EO +év)-F(y)| 
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where, s* is the increment of the kth Newton iteration. 
By replacing |F "y| and y in Eq. (9) with that given in Eq. (10) and (11), the optimal perturbation size 


is computed using Eq. (12). 
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From the above derivation, it can be observed that there is almost no additional computing cost in the estimation 


of the optimal perturbation size. Additionally, the ratio of M / y is the most important, not the value of M or 
y itself. 


3 Implementation of coupling neutronics with thermal-hydraulics 


3.1 High-fidelity Neutron Transport Program (HNET) 

HNET is a 3D whole-core transport code capable of performing high-fidelity neutron transport calculations for 
both steady-state and transient problems. The 2D/1D method with multilevel generalized equivalence theory-based 
CMEFED acceleration (gCMFD) was developed in HNET for transport calculations. Specifically, the 2D method of 
characteristics (MOC) was utilized for the radial solver, and the two-node nodal expansion method is used for the 
axial dependency. The 2D and 1D solutions were coupled using transverse leakage. To decrease the computational 
cost of the CMFD, an innovative self-developed RSILU-preconditioned GMRES solver was utilized for the CMFD 
linear system. Regarding transient simulation, the multi-level predictor-corrector quasi-static scheme, together with 
the multi-level gCMFD acceleration, was developed to accelerate the time-dependent neutron transport calculation. 
Additionally, HNET employs the 47-group HELIOS library to provide cross sections, and the subgroup method was 
used for resonance self-shielding calculations. A detailed description of the transport calculations in HNET can be 
found in the References [14, 15, 19, 20]. 

3.2 Thermal-hydraulic model 

To provide the fuel temperature, coolant temperature, and coolant density to the neutronics field, a simplified 
internal thermal-hydraulic module was developed in the HNET. For pressurized water reactors, the module consists 
of a 1D radial heat conduction model and a 1D axial heat convection model, which have been proven to be effective 
for a wide range of applications [8]. 

The heat-transfer calculation began with a cylindrical conduction equation. Because of the large ratio of axil 
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length to radius and the axisymmetric geometry, axial and circumferential heat transfer can be ignored; thus, only 
radial heat conduction is considered for the fuel rods, as in Eq. (13). 


oT 10 oT 
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where T, p, Cp, k, and r denote temperature, density, specific heat, thermal conductivity, and radius, respectively. The 
volumetric power density q” contributes to the fission reaction. 


In addition, it is necessary to provide boundary conditions between heat conduction and convection. The heat 
flux transferred from the pellet to the fluid is determined by the following boundary condition as in Eq. (14). 


a es eee (14) 
i or 
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where Tm and 7; are the temperatures of the cladding outer wall and coolant, respectively, and hy, is the convective 
heat transfer coefficient, which is determined by the Dittus-Boelter correlation. 

In terms of the flow field, a single-channel model with the following reasonable assumptions is applied: (a) 
Each assembly is defined as a channel region. In addition, no cross-flow between the adjacent channels was 
considered. (b) The coolant was always in the liquid phase. (c) The coolant pressure was approximately constant 
during the steady state and was fast transient. 

Under these assumptions, the momentum conservation equation was removed. Thus, only the remaining mass 
and energy conservation equations were resolved in the flow field, and their general formulas in the axial direction 
are given by Eq. (15) and Eq. (16). 
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where p, v, and h denote the density, velocity, and enthalpy of coolant, respectively. z is the height in the axial 
direction and ¢ represents time. Ly and A; are the heated perimeter and the cross-sectional area of the channel, 


respectively. qi. is the heat flux at the cladding surface, and q” is the heat source generated directly inside the 


coolant. 

To solve the 1D radial heat conduction equations and 1D axial heat convection equations, the finite volume 
method is used to carry out spatial discretization, the theta method is adopted to discretize the time derivative terms, 
and the discrete formulas for the numerical calculation can be obtained. 

3.3 Coupling implementation 

The implementation of the MFNK method in the steady-state coupling calculation is shown in Fig. 1. Notably, 
two coupling configurations can be arranged according to the choice of the initial guess. Because the dimension of 
power density is small in comparison with thermal-hydraulic solutions, the power density was selected as the initial 
estimate in this study. Next, to initiate the MFNK scheme with an appropriate initial power density y, the Picard 
iteration was performed first. Specifically, according to a primary y, thermal-hydraulic solutions are calculated and 


then used to update the macroscopic cross sections. Then, based on the results of the transport calculation, a new 
Picard 


y was obtained. 
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Compute Ay’ = Ay? +Vq; 


Fig. 1 Solution scheme of MFNK for steady-state coupling calculation 


The MFNK method began after the initial guess calculation. Based on the given y, the Picard iteration was 
performed to form the residual, as described in Eq. (3). Subsequently, the Newton iteration was solved using the 
GMRES method, where the required matrix-vector product was approximated through a finite difference 
approximation with the optimal perturbation size. In each Newton step, a new solution y was updated. The iterative 
scheme was repeated until convergence was achieved. In addition, the maximum Krylov subspace dimension was 
chosen as two, and the transport solver was fully converged in each calculation. 

In the transient simulation, the OS method was used to implement the transient coupling scheme. To improve 
stability, a staggered time mesh was applied to the OS method [21]. In this scheme, the two physical fields are 
advanced for different time step sizes. The thermal-hydraulics problem was solved first with one-half of the 
neutronics time step. Subsequently, based on the updated thermal-hydraulic conditions, the neutronics solver was 
invoked to calculate the neutron flux. Then, the thermal-hydraulics solver was advanced again using the new power 


density with one-half of the neutronics time step. 
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Fig. 2 (Color online) OS method with staggered time step 


4 Numerical results and analysis 


VERA Problem 6 [17] with different power levels and boron concentrations was first employed to study the 
convergence performance of the MFNK method. A comparison of the numerical results for MFNK, Picard, and 
Picard with relaxation is presented in this section. Subsequently, the NEACRP transient benchmark [18] was used 
to preliminarily test and verify the transient coupling calculation. The corresponding results were compared with 
reference solutions for verification. 

4.1 Assessment of MFNK performance 

VERA Problem 6 describes a single fuel assembly problem under typical full-power and nominal thermal- 
hydraulic conditions, which are presented in Table 1. The assembly consisted of 264 fuel rods, 24 guide tubes, and 
an instrument tube at the center. Meanwhile, there were six intermediate grids and two end grids in the assembly 
along the axial direction. More detailed geometric and material specifications are provided in [17]. 

Table 1 Nominal thermal-hydraulics conditions 


Parameter Value 
System pressure (psia) 2250 
Rated power (100%) (MW) 17.67 
Rated coolant mass flow (Mlbs/hr) 0.6823 
Inlet coolant temperature (K) 565 
Boron concentration (ppm) 1300 


For VERA Problem 6, the transport equation was solved using the 2D/1D method via gCMFD acceleration. 
Based on the 47-group HELIOS library, nuclide density, and temperature, the cross-sections were obtained through 
resonance treatment. Fig. 3 illustrates the discretization of MOC mesh in radial direction, the heterogeneous fuel 
pins are divided into 56 flat source regions (FSR), and the reflector cells use the 2 x 2 Cartesian sub-mesh. In 
addition, a ray spacing of 0.03 cm with 64 azimuthal angles and a Tabuchi- Yamamoto polar quadrature using 3 polar 
angles per half-space were utilized for MOC calculation. The active fuel region was divided into 49 axial layers, 
similar to VERA Problem 3 [17]. Furthermore, the core plates, nozzles, and grids were explicitly modeled for the 
transport domain. The radial boundaries were set to be reflective, and the vacuum boundary condition was applied 
axially, owing to the top and bottom water reflectors. For the thermal-hydraulic domain, each fuel pin, including the 
fuel pellet, gap, and cladding, was explicitly modeled and calculated individually using the radial power density 
inside the fuel pin. The discretization of the coolant channel in the axial direction is the same as that in the transport 
calculation. In addition, the inlet coolant conditions were adopted for the bottom structural material below the active 
fuel, whereas the average exit conditions of the coolant at the top of the active fuel region were used as the top 


structural material. 
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(a) (b) 
Fig. 3 (Color online) The division of MOC mesh in radial direction: (a) fuel pin; (b) reflector cell 
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Generally, a fission source is used to measure the convergence of a coupled system. To evaluate the 
convergence behavior of MFNK, a strict convergence criterion using the power density residual norm was set for 
VERA Problem 6, as in Eq. (17). 


2 
RMS = FOV 21x10% (17) 
n 


where RMS is the root mean square of the residual vector and n is the linear system dimension. 

First, according to Eq. (6), the influence of perturbation size on the MFNK method was investigated through 
VERA Problem 6. To test the behavior of the optimal perturbation size defined by Eq. (12), three empirical 
perturbation sizes were considered for the comparison. The first was a simple empirical perturbation with a constant 
value of 0.001. The second method is based on Pernice’s method [22]. 
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The third empirical perturbation size proposed by Knoll [11] is given by Eq. (19). 
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where b is a constant whose magnitude is within a few orders of magnitude of the square root of the machine round- 
off. 

Fig. 4 illustrates how the perturbation size affects the convergence behavior of the MFNK method for VERA 
Problem 6. It can be observed that MFNK with a constant perturbation size converges very slowly. The 
performances of the other two empirical perturbation sizes are relatively consistent, and their residual norms seem 
to stagnate around 0.1. By contrast, MFNK with an optimal perturbation size can converge further and satisfy the 
convergence criterion. The primary reason is that the optimal perturbation size can reduce the error introduced by 
the finite-difference approximation even when the residual norm becomes small, whereas the empirical ones cannot. 


Hence, the optimal perturbation size is used in the following analysis. 


Power density residual norm 


1 
l 
1 —=— Optimal 
10° —*— Pernice 
—— Knoll 
as —— 0.001 
0 5 10 15 20 


MENK iteration index 
Fig. 4 The effect of perturbation size in the MFNK method 


Subsequently, a convergence analysis was performed by comparing the results between the MFNK and Picard 
methods. Notably, in order to make an appropriate comparison between the MFNK and Picard methods, the 
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convergence history was monitored by the iteration index of the MOC, which is a time-consuming procedure, rather 
than their individual outer iterations. 

The convergence history of ker and the maximum fuel temperature for VERA problem 6 are shown in Fig. 5 
and Fig. 6, respectively. The results indicate that MFNK achieved a stable state quickly, whereas Picard initially 
exhibited oscillatory behavior. For instance, the oscillation of keg in MFNK was a maximum of 1 pem, whereas the 
maximum variation of kerr in Picard is approximately 38 pcm. In addition, as illustrated in Fig. 5 and Fig. 6, the 
solutions resulting from MFNK and Picard match well. 
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Fig. 5 (Color online) Convergence history for ker 
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Fig. 6 (Color online) Convergence history for maximum fuel temperature 


To further evaluate the performance of MFNK, VERA Problem 6 with a modification of the power level and 
boron concentration was simulated. Traditionally, relaxation has been introduced to Picard for better convergence 
behavior. Hence, Picard and Picard with relaxation applied to the temperature solution were performed for 
comparison. Because the optimal choice of relaxation factor depends strongly on the problems, a relaxation factor 
of w = 0.5 was adopted in this paper [8]. The overall convergence performances of MFNK, Picard, and Picard with 
relaxation are summarized in Table 2 and Table 3. In addition to the MOC iterations required for convergence, the 
total number of resonance calculations is also presented because of the relatively high computational cost. Notably, 
the resonance calculation was performed in the coupling scheme only when the maximum temperature change 
between two adjacent resonance calculations was larger than 1 K. 

As shown in Table 2, MFNK and Picard with relaxation can considerably reduce the MOC iterations compared 
with Picard over all except for the 60% power level case. Although the three methods displayed a dependence on 
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the power magnitude, MFNK and Picard with relaxation were less susceptible to the power level than Picard. 
Furthermore, it was observed that MFNK and Picard with relaxation could achieve convergence criteria in all cases, 
whereas Picard performed poorly and even diverged as the power level increased. This is primarily because 
increasing the power enhances coupling strength. However, Picard is generally unable to handle this stronger 
coupling. Additionally, the number of MOC iterations for MFNK was the same as or less than that for Picard with 
relaxation. More specifically, MFNK can decrease the number of MOC iterations from Picard with a maximum 
relaxation of 15. 

In terms of boron perturbation, Table 3 demonstrates that varying the boron concentration does not affect the 
convergence behavior. As shown in Table 3, the performance of MFNK is slightly better than that of Picard with 
relaxation, and poor convergence can be found in the Picard results. Table 2 and Table 3 show that the number of 
resonance calculations for MFNK and Picard with relaxation are comparable, and both are less than that for Picard 
as a whole. In other words, the cross sections can converge faster in MFNK and Picard with relaxation. 

Table 2 Iterations to convergence for VERA Problem 6 at various power levels 


Power level (%) 


Couple Scheme Iterations 
60 80 100 120 140 
MOC 85 105 107 126 122 
MFNK 
Resonance 10 12 11 10 11 
; MOC 103 113 117 126 137 
Picard (w = 0.5) 
Resonance 11 11 11 12 12 
Picard MOC 72 122 197 divergence divergence 
icar 
Resonance 5 10 26 divergence divergence 


Table 3 Iterations to convergence for VERA Problem 6 at various boron concentrations 


Boron concentration (ppm) 


Couple Scheme Iterations 
0 500 900 1300 1700 
MOC 112 115 111 107 113 
MFNK 
Resonance 11 12 12 11 13 
A MOC 119 116 117 117 122 
Picard (w = 0.5) 
Resonance 12 12 12 11 11 
. MOC 177 209 223 197 216 
Picard 
Resonance 20 27 30 26 29 


Fig. 7 compares the power density residual norm of MFNK against Picard with and without relaxation with 
respect to the variation in the power level. In contrast to Picard with relaxation, the residual norm of MFNK dropped 
more rapidly in the 80% and 140% power cases. In addition, Picard converged at a markedly slower rate and 
divergence behavior occurred in Picard at a 140% power level. 


—=— MFNK 


—— Picard 
—— Picard (w = 0.5) 


—— MFNK 
—— Picard 
—— Picard (w = 0.5) 


Power density residual norm 
Power density residual norm 


l 1 
1 l 
1 1 
1 1 
1 10 
106 10° 
40 60 80 100 120 40 60 80 100 120 140 160 180 200 
MOC yee index MOC Se index 
a 


Fig. 7 (Color online) Convergence behavior of the residual norm for VERA Problem 6 with power perturbation: 
(a) 80% power level; (b) 140% power level 


Based on the above analysis, it can be concluded that although MFNK does not perform as efficiently as Picard 
at the 60% power level, it can display a better convergence performance than Picard for a wider range of problems. 
Moreover, although the performance of MFNK is better than that of Picard with a relaxation factor of 0.5 in given 
cases, the opposite results may occur if the optimum values of the relaxation factor are selected. Nevertheless, it 
should be easier to improve the convergence rate and enhance the numerical stability when utilizing MFNK than 
when applying a relaxation scheme, because the optimal relaxation factor is problem-dependent and generally 
chosen empirically [21]. 


4.2 Transient simulation of control rod ejection 

The NEACRP benchmark is a 3D PWR core transient problem constructed to simulate the ejection of a control 
assembly (CA) from the core under operational conditions. The assembly arrangement is shown in Fig. 8. Moreover, 
2-group macroscopic cross sections were set for 11 homogenized compositions. However, the NEACRP benchmark 
was originally specified for diffusion-based codes, resulting in a lack of total and self-scattering cross sections. To 
obtain an equivalent transport version, the total cross sections can be approximated to transport cross sections, and 


then the self-scattering cross sections are calculated using Eq. (20). 


G 

Eog “Ere Eag- 2 Poy (20) 
oe 

where Ł, Ya, and È, are the total, absorption, and scattering cross-sections, respectively. 

The NEACRP benchmark comprises six cases with different combinations of control arrangements, original 
positions, and power levels. A detailed description of the geometry, macroscopic cross-sections, and initial 
configuration of the control assemblies can be found in [18]. In this study, cases A2 and B2 were considered to 
verify the accuracy of the transient simulation. Both cases were initiated by the rapid ejection of the CA in 0.1 s at 
hot full power state; the CA to be ejected is marked by a circle in Fig. 8. Specifically, the symbols A and B represent 
the ejected CA of cases A2 and B2, respectively. 
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Fig. 8 (Color online) Radial layout of the quadrant core in NEACRP benchmark 
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To resolve the NEACRP benchmark using the transient 2D/1D transport solver in HNET [15], the geometry 
was uniformly discretized into an FSR mesh size of 0.7202 cm x 0.7202 cm for the MOC calculation. Additionally, 
the ray spacing for the NEACRP benchmark is the same as that for VERA Problem 6. From bottom to top, the active 
zone was divided into 26 layers with heights of 7.7, 11, 15 (21 layers), 12.8 (2 layers), and 8.0 cm, which differed 
from the division of the NEACRP benchmark for high accuracy of results. A reflective boundary was used to the 
west and north because of symmetry, whereas the rest used vacuum boundary conditions. Concerning the thermal- 
hydraulic calculation, an average fuel pin, representing all fuel pins of the associated assembly, was modeled in 
detail to obtain the radial temperature distribution by applying the average assembly power density. In the axial 
direction, the separation of the coolant channel was consistent with that of the transport calculation. 

The steady-state thermal-hydraulic conditions are listed in Table 4. The relation between thermal properties, 
such as thermal conductivity and specific heat capacity, and temperature is derived from [18]. Moreover, to consider 
the feedback effects from thermal hydraulics and boron concentration, the macroscopic cross-sections for each 
composition are defined by Eq. (21). 


se 5 (Sea) (21) 
= Pn IF Cn 
where p, Tm, Te and C are the moderator temperature, moderator density, fuel temperature, and boron density, 
respectively. The subscript 0 represents the reference values. 
Table 4 Steady-state thermal-hydraulic characteristics 


Parameters Value 

Core thermal output (MW) 2775 (HFP) 
Core pressure (Mpa) 15.5 
Coolant inlet temperature (°C) 286 

Mass flow (kg/s) 12893 


A comparison of steady-state and transient solutions for cases A2 and B2 is given in Table 5. The reference 
results listed in Table 5 were derived from [23]. As shown in Table 5, there is good agreement between the simulation 
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and reference results. The critical boron concentration at the steady state deviated slightly from the reference value, 


and the maximum difference was approximately 31 ppm, which is considered acceptable. In addition, the 


representative parameters in the transient process coincided well with the reference. 
Table 5 Steady-state and transient calculation results 


A2 B2 
Parameters 


Reference Calculation Deviation Reference Calculation Deviation 


Critical boron 


: 1156.63 1187.66 31.03 1183.83 1214.79 
concentration (ppm) 
Power peak 1.0800 1.0821 0.0021 1.0630 1.0651 
Peak time (s) 0.0100 0.0105 0.0005 0.1200 0.1200 
Final power at 5 s 1.0350 1.0357 0.0007 1.0380 1.0399 
Coolant outlet average 
324.600 325.049 0.449 324.700 325.164 


temperature at 5 s (°C) 


Fig. 9 depicts the evolution of the core-normalized power and reactivity for case A2. Here, 


30.96 


0.0021 
0 
0.0019 


0.464 


the adopted 


reference solutions are computed using the PARCS code [24]. As indicated, the HNET results for case A2 are in 
agreement with that of PARCS. The slight deviation in the normalized power and reactivity curves can be explained 
by the difference in the neutronics methodology. The NEACRP benchmark is resolved by neutron transport 


calculations in HNET, whereas the nodal method is used in PARCS. In addition, as illustrated in Fig. 


10 (a), there 


is no significant difference in the responses of the coolant exit temperature to the CA ejection between these two 
codes. Because the result of the Doppler temperature from the PARCS is not provided in [24], Fig. 10 (b) plots only 
the solution of the HNET. From Fig. 10 (b), an increase in the average Doppler temperature owing to power 


increments is observed, which behaves as expected. 
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Fig. 9 (Color online) Transient responses in case A2: (a) normalized power; (b) reactivity 
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Fig. 10 (Color online) Transient responses in case A2: (a) coolant exit temperature; (b) core averaged Doppler 


temperature 


As shown in Fig. 11 and Fig. 12, to further test the level of agreement, a closer analysis of the axial relative 
power and radial normalized power for case A2 was performed at 0 s, 0.1 s, and 5 s. The relative error e is defined 
by Eq. (22). 

calc — re 
e= One id x 100% (22) 
ref 
where calc represents the values calculated by HNET, and ref'is the reference value derived from PARCS [24]. 

The normalized power in the pin and assembly levels are calculated using Eqs. (23) and (24). 


= PN. 

P pin = ET 7 
D P oon N 

P assem = _assem___assem D 


core 


where Ppin and Passem denote the axial-integrated power per fuel pin and per fuel assembly, respectively; Npin and 
Nassem represent the total number of fuel pins and fuel assemblies, respectively; and Peore is the core power. 

With respect to the axial relative power, no significant difference between HNET and PARCS is observed in 
Fig. 11. Fig. 12 shows that the results of assembly level normalized power at the considered temporal points from 
HNET agree well with PARCS. Moreover, because the NEACRP benchmark was performed using a transient 
transport solver, the pin-wise power distribution can be obtained directly without power reconstruction. Fig. 13 
depicts the pin-wise radial normalized power for case A2 at 0.0 s and 0.1 s. The comparison shows that the pin 
power in the center increases significantly after the withdrawal of the control rod, which coincides with the sharp 
increase in power observed in Fig. 9. The above results prove that the evolution of the power distribution can be 
effectively captured when the reactivity is being inserted. 
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Fig. 11 (Color online) Axial relative power for case A2: (a) time at 0.0 s; (b) time at 0.1 s; (c) time at 5.0 s 
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Fig. 12 (Color online) Relative error of assembly normalized power for case A2 in quarter core geometry (%): (a) 


time at 0.0 s; (b) time at 0.1 s; (c) time at 5.0 s 


(a) (b) 


Fig. 13 (Color online) Pin normalized power distribution for case A2 in quarter core geometry: (a) time at 0.0 s; 
(b) time at 0.1 s 
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Finally, because the thermal-hydraulics feedback is provided at the assembly level for the NEACRP benchmark, 
Fig. 14 shows the local temperatures in assemblies C-3 and B-8, where the highest and lowest powers appear, 
respectively. It can be observed that there are obvious discrepancies for fuel centerline temperatures between the 
two assemblies. As shown in Fig. 14 (b), the increase in temperatures between 0 s and 5 s in assembly B-8 is larger 
than that in assembly C-3, particularly the fuel centerline temperature. 
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Fig. 14 (Color online) Local temperatures and change in temperatures in the assemblies with the highest and 


lowest powers: (a) local temperatures; (b) change in temperatures between 0 s and 5 s 
5 Conclusion 


The primary objective of this study was to couple neutronics with thermal hydraulics in HNET. First, a 
simplified internal thermal-hydraulic module was developed for thermal properties. Subsequently, a steady-state 
coupling calculation is achieved by implementing the MFNK method. Next, the OS method with a staggered time 
mesh was used to accomplish the transient coupling simulation. 

The investigation of the perturbation size shows that the optimal perturbation size can significantly enhance 
the convergence of the MFNK method compared with the empirical perturbation size. With respect to VERA 
problem 6 at different power levels and boron concentrations, the MFNK method is slightly sensitive to the feedback 
intensity and can considerably reduce the computational cost compared with the Picard iteration in most cases. It 
can be proven that the MFNK method is more stable and efficient than the Picard iteration for a wider range of 
problems. Notably, the performance of MFNK is as good as or better than that of Picard with relaxation, and there 
is no requirement for input data in MFNK compared with the relaxation approach. Furthermore, the numerical 
results for the NEACRP transient benchmark are in agreement with the reference values, indicating the validation 
of the transient coupling calculation in the hot full-power state. In conclusion, the results demonstrate that the 


neutronics and thermal-hydraulics coupling simulation is effectively performed in HNET under both steady-state 
and transient conditions. 
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